#script to run model for photosynthesis data – similar to growth_rate script
library(lme4)
## Loading required package: Matrix
library(lmerTest)
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(afex)
## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: emmeans() for ANOVA models now uses model = 'multivariate' as default.
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
##
## Attaching package: 'afex'
## The following object is masked from 'package:lme4':
##
## lmer
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(car)
library(MuMIn)
library (dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
##
## recode
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(emmeans)
library(DHARMa)
## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(performance)
library(patchwork)
library(rstatix)
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2
## ──
## ✔ tibble 3.1.8 ✔ purrr 1.0.1
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.3
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ rstatix::filter() masks dplyr::filter(), stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some() masks car::some()
## ✖ tidyr::unpack() masks Matrix::unpack()
library(sjPlot)
## Learn more about sjPlot with 'browseVignettes("sjPlot")'.
library(sjmisc)
##
## Attaching package: 'sjmisc'
##
## The following object is masked from 'package:purrr':
##
## is_empty
##
## The following object is masked from 'package:tidyr':
##
## replace_na
##
## The following object is masked from 'package:tibble':
##
## add_case
library(mmtable2)
##
## Attaching package: 'mmtable2'
##
## The following object is masked from 'package:tidyr':
##
## table1
library(gt)
library(purrr)
library(stringr)
library(tidyr)
#load this file for rETR-based use per Silsbe and Kromkamp 2012 #all_runs_photosyn_data <- read.csv(“data_input/hyp_ulva_all_runs_ek_alpha.csv”)
#load this file for normalized to quantum efficiency of photosynthesis per same as above
all_runs_photosyn_data <- read.csv("data_input/hyp_ulva_all_runs_ek_alpha_normalized.csv")
all_runs_photosyn_data$Run <- as.factor(all_runs_photosyn_data$Run)
#assign temperature as a factor
all_runs_photosyn_data$Temperature <- as.factor(all_runs_photosyn_data$Temp...C.)
#assigns treatment as characters from integers then to factors
all_runs_photosyn_data$Treatment <- as.factor(as.character(all_runs_photosyn_data$Treatment))
all_runs_photosyn_data$deltaNPQ <- as.factor(all_runs_photosyn_data$deltaNPQ)
#toggle between the species for output. Use Day 9 for final analysis
ulva <- subset(all_runs_photosyn_data, Species == "ul" & RLC.Day == 9 & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol"
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol"
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol"
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"
#ULVA pmax________________________________________________________________
ulva %>% ggplot(aes(pmax)) +
geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
#run model without interaction between the treatments and temperature
ulva_pmax_model <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)
#construct null model to perform likelihood ratio test REML must be FALSE
ulva_pmax_treatment_null <- lmer(formula = pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
## boundary (singular) fit: see help('isSingular')
ulva_pmax_model2 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_pmax_treatment_null, ulva_pmax_model2)
## Data: ulva
## Models:
## ulva_pmax_treatment_null: pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_pmax_model2: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## ulva_pmax_treatment_null 7 2042.7 2067.1 -1014.37 2028.7
## ulva_pmax_model2 11 1919.9 1958.2 -948.94 1897.9 130.86 4
## Pr(>Chisq)
## ulva_pmax_treatment_null
## ulva_pmax_model2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ulva_pmax_temperature_null <- lmer(formula = pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_pmax_model3 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_pmax_temperature_null, ulva_pmax_model3)
## Data: ulva
## Models:
## ulva_pmax_temperature_null: pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_pmax_model3: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## ulva_pmax_temperature_null 9 1919.7 1951.0 -950.83 1901.7
## ulva_pmax_model3 11 1919.9 1958.2 -948.94 1897.9 3.778 2
## Pr(>Chisq)
## ulva_pmax_temperature_null
## ulva_pmax_model3 0.1512
#make residual plots of the data for ulva
hist(resid(ulva_pmax_model))
plot(resid(ulva_pmax_model) ~ fitted(ulva_pmax_model))
qqnorm(resid(ulva_pmax_model))
qqline(resid(ulva_pmax_model))
#check the performance of the model
performance ::check_model(ulva_pmax_model)
## Variable `Component` is not in your data frame :/
r.squaredGLMM(ulva_pmax_model)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
## R2m R2c
## [1,] 0.5918687 0.6885467
tab_model(ulva_pmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
| pmax | ||||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | 34.05 | 3.31 | 27.52 – 40.57 | 10.28 | <0.001 | 229.00 |
| Treatment [1] | 21.54 | 3.54 | 14.57 – 28.51 | 6.09 | <0.001 | 229.00 |
| Treatment [2] | 22.52 | 3.54 | 15.55 – 29.49 | 6.36 | <0.001 | 229.00 |
| Treatment [3] | 41.51 | 3.54 | 34.54 – 48.49 | 11.73 | <0.001 | 229.00 |
| Treatment [4] | 43.48 | 3.54 | 36.50 – 50.45 | 12.28 | <0.001 | 229.00 |
| Temperature [27] | -5.55 | 2.80 | -11.06 – -0.04 | -1.99 | 0.048 | 229.00 |
| Temperature [30] | -3.22 | 2.54 | -8.22 – 1.78 | -1.27 | 0.205 | 229.00 |
| Random Effects | ||||||
| σ2 | 135.25 | |||||
| τ00 Plant.ID | 26.61 | |||||
| τ00 Run | 8.13 | |||||
| τ00 RLC.Order | 7.25 | |||||
| ICC | 0.24 | |||||
| N Run | 7 | |||||
| N Plant.ID | 95 | |||||
| N RLC.Order | 6 | |||||
| Observations | 240 | |||||
| Marginal R2 / Conditional R2 | 0.592 / 0.689 | |||||
plot(allEffects(ulva_pmax_model))
ulva %>% ggplot(aes(treatment_graph, pmax)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) +
labs(x="salinity/nitrate", y= "Day 9 Pmax (μmols electrons m-2 s-1)", title= "A", subtitle = "Ulva lactuca") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(-1, 170) + stat_mean() +
scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
geom_hline(yintercept=50, color = "red", size = 0.5, alpha = 0.5) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
#summarize the means for pmax
ulva %>% group_by(Treatment) %>% summarise_at(vars(pmax), list(mean = mean))
## # A tibble: 5 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 31.1
## 2 1 52.3
## 3 2 53.3
## 4 3 72.3
## 5 4 74.3
ulva %>% group_by(Run) %>% summarise_at(vars(pmax), list(mean = mean))
## # A tibble: 7 × 2
## Run mean
## <fct> <dbl>
## 1 1 66.3
## 2 2 58.9
## 3 3 64.7
## 4 3.5 67.5
## 5 4 60.9
## 6 6 34.2
## 7 6.5 28.1
ulva %>% group_by(Treatment, RLC.Day) %>% summarise_at(vars(pmax), list(mean = mean))
## # A tibble: 5 × 3
## # Groups: Treatment [5]
## Treatment RLC.Day mean
## <fct> <int> <dbl>
## 1 0 9 31.1
## 2 1 9 52.3
## 3 2 9 53.3
## 4 3 9 72.3
## 5 4 9 74.3
#add growth rate from other dataset to this one and subset by species
growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$lunar.phase <- as.factor(growth_rate$lunar.phase)
#make a new column for weight change (difference final from initial)
growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100
gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)
ulva$lunar.phase <- (gr_ulva$lunar.phase)
#plot a regression between the photosynthetic independent variables of interest and growth rate
ulva_growth_etr_graph <- ggplot(ulva, aes(x=pmax, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "Ulva lactuca Pmax vs Growth Rate", x = "Pmax (μmols electrons m-2 s-1)",
y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_etr_graph
## `geom_smooth()` using formula = 'y ~ x'
#_____________________________________________________________________________________________________________ #HYPNEA
#There was no D9 RLC for hm6-4 on 11/12/21 but had to remove hm6-4 from 10/9/21 below to match growth data
hypnea <- subset(all_runs_photosyn_data, Species == "hm" & RLC.Day == 9 & uid != "2021-10-09_hm6-4")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol"
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol"
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol"
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"
#for Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) # and hm6-4 on 10/9/21 because it was white and also looked dead
gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)
#hm1-1 on 10/12/22 and hm1-2 on 4/29/22 causing issues of influential observations #hm1-2 for both pmax and Ek – leaving them in dataset because no good reason to believe not good data #make a histogram and residual plots of the data for hypnea
hist(hypnea$pmax, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)
hypnea %>% ggplot(aes(pmax)) +
geom_histogram(binwidth=5, fill = "#7D0033", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
#run model for pmax
hyp_pmax_model <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea)
hist(resid(hyp_pmax_model))
plot(resid(hyp_pmax_model) ~ fitted(hyp_pmax_model))
qqnorm(resid(hyp_pmax_model))
qqline(resid(hyp_pmax_model))
#check the performance of the model
performance::check_model(hyp_pmax_model)
## Variable `Component` is not in your data frame :/
r.squaredGLMM(hyp_pmax_model)
## R2m R2c
## [1,] 0.2732395 0.6374791
summary(hyp_pmax_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
## (1 | RLC.Order)
## Data: hypnea
##
## REML criterion at convergence: 2271.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4349 -0.4862 -0.0761 0.4488 3.7551
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 42.71 6.535
## Run (Intercept) 83.40 9.132
## RLC.Order (Intercept) 14.42 3.797
## Residual 139.86 11.826
## Number of obs: 286, groups: Plant.ID, 143; Run, 8; RLC.Order, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 56.175 7.120 6.156 7.890 0.000193 ***
## Treatment1 2.421 8.135 5.388 0.298 0.777153
## Treatment2 -1.295 8.135 5.388 -0.159 0.879349
## Treatment2.5 14.733 11.520 4.722 1.279 0.260159
## Treatment3 3.905 8.135 5.388 0.480 0.650025
## Treatment4 7.689 8.143 5.411 0.944 0.385299
## Temperature27 -17.253 3.116 32.273 -5.537 4.05e-06 ***
## Temperature30 -19.429 2.729 52.871 -7.119 2.94e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.786
## Treatment2 -0.786 0.956
## Treatmnt2.5 -0.555 0.486 0.486
## Treatment3 -0.786 0.956 0.956 0.486
## Treatment4 -0.785 0.955 0.955 0.485 0.955
## Temperatr27 -0.206 0.001 0.001 0.000 0.001 0.000
## Temperatr30 -0.196 0.000 0.000 0.000 0.000 0.003 0.467
tab_model(hyp_pmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
| pmax | ||||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | 56.18 | 7.12 | 42.16 – 70.19 | 7.89 | <0.001 | 274.00 |
| Treatment [1] | 2.42 | 8.13 | -13.59 – 18.44 | 0.30 | 0.766 | 274.00 |
| Treatment [2] | -1.29 | 8.13 | -17.31 – 14.72 | -0.16 | 0.874 | 274.00 |
| Treatment [2.5] | 14.73 | 11.52 | -7.95 – 37.41 | 1.28 | 0.202 | 274.00 |
| Treatment [3] | 3.91 | 8.13 | -12.11 – 19.92 | 0.48 | 0.632 | 274.00 |
| Treatment [4] | 7.69 | 8.14 | -8.34 – 23.72 | 0.94 | 0.346 | 274.00 |
| Temperature [27] | -17.25 | 3.12 | -23.39 – -11.12 | -5.54 | <0.001 | 274.00 |
| Temperature [30] | -19.43 | 2.73 | -24.80 – -14.06 | -7.12 | <0.001 | 274.00 |
| Random Effects | ||||||
| σ2 | 139.86 | |||||
| τ00 Plant.ID | 42.71 | |||||
| τ00 Run | 83.40 | |||||
| τ00 RLC.Order | 14.42 | |||||
| ICC | 0.50 | |||||
| N Run | 8 | |||||
| N Plant.ID | 143 | |||||
| N RLC.Order | 6 | |||||
| Observations | 286 | |||||
| Marginal R2 / Conditional R2 | 0.273 / 0.637 | |||||
plot(allEffects(hyp_pmax_model))
#construct null model to perform likelihood ratio test REML must be FALSE
hypnea_pmax_treatment_null <- lmer(formula = pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_pmax_model2 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_pmax_treatment_null, hypnea_pmax_model2)
## Data: hypnea
## Models:
## hypnea_pmax_treatment_null: pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_pmax_model2: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## hypnea_pmax_treatment_null 7 2335.5 2361.1 -1160.8 2321.5
## hypnea_pmax_model2 12 2329.8 2373.7 -1152.9 2305.8 15.755 5
## Pr(>Chisq)
## hypnea_pmax_treatment_null
## hypnea_pmax_model2 0.007581 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hypnea_pmax_temperature_null <- lmer(formula = pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_pmax_model3 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_pmax_temperature_null, hypnea_pmax_model3)
## Data: hypnea
## Models:
## hypnea_pmax_temperature_null: pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_pmax_model3: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## hypnea_pmax_temperature_null 10 2365.9 2402.4 -1173.0 2345.9
## hypnea_pmax_model3 12 2329.8 2373.7 -1152.9 2305.8 40.109 2
## Pr(>Chisq)
## hypnea_pmax_temperature_null
## hypnea_pmax_model3 1.952e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hypnea %>% ggplot(aes(treatment_graph, pmax)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) +
labs(x="Treatment", y= "Day 9 Pmax (μmols electrons m-2 s-1)", title= "B", subtitle = "Hypnea musciformis") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(-1, 170) + stat_mean() +
scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
geom_hline(yintercept=50, color = "red", size = 0.5, alpha = 0.5) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
#plot a regression between the photosynthetic independent variables of interest and growth rate
hypnea_growth_etr_graph <- ggplot(hypnea, aes(x=pmax, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "Hypnea musciformis Pmax vs Growth Rate", x = "Pmax (μmols electrons m-2 s-1)",
y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor(label.y = 150)
hypnea_growth_etr_graph
## `geom_smooth()` using formula = 'y ~ x'
#summarize the means
hypnea %>% group_by(Treatment) %>% summarise_at(vars(pmax), list(mean = mean))
## # A tibble: 6 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 43.9
## 2 1 47.5
## 3 2 43.8
## 4 2.5 58.7
## 5 3 49.0
## 6 4 53.1
hypnea %>% group_by(Treatment, RLC.Day) %>% summarise_at(vars(pmax), list(mean = mean))
## # A tibble: 6 × 3
## # Groups: Treatment [6]
## Treatment RLC.Day mean
## <fct> <int> <dbl>
## 1 0 9 43.9
## 2 1 9 47.5
## 3 2 9 43.8
## 4 2.5 9 58.7
## 5 3 9 49.0
## 6 4 9 53.1
#plot a regression between pmax and ek
hypnea_etr_ek_graph <- ggplot(hypnea, aes(x=pmax, y=ek.1)) +
geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "Hypnea musciformis pmax vs Ek", x = "Pmax (μmols electrons m-2 s-1)",
y = "Ek (μmols photons m-2 s-1)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
hypnea_etr_ek_graph
## `geom_smooth()` using formula = 'y ~ x'